% Function : Forward Search Algorithm with Correlation(extended version)
% Parameter:
%   Y : target data that we are trying to predict
%   X : training set, each row is one training example, each column is one feature
%   u : upper limit of selected features, search depth
%   gc: values of gray correlation
%   s : separator, [1...s] are training examples, [s+1...end] are validation ones
% Return:
%   I: indexes of selected features
%   R : average residuals of selected features
%   RR: R^2 statistic of selected features
function [I, R, RR] = CForwardx(Y, X, u, gc, s)
    if nargin ~= 5
        error('Usage: Forward(Y, X, u, gc, s)');
    end

    [m, ~] = size(X);   % m: number of samples, n: number of features

    [~, L] = KMeans(gc, 2);

    % initialize collection
    F  = ones(m, u + 1);    % collection of features
    I  = zeros(1, u);       % index of collected features
    R  = zeros(1, u);       % average residuals
    RR = zeros(1, u);       % R^2 statistic

    i = 0;
    k = 1;
    while k <= u
        i = mod(i, 3) + 1;  % i'th cluster

        % indexes of features in i'th cluster
        ft  = find(L == i);
        nft = numel(ft);    % number of features

        % r  = zeros(1, nft);
        % rr = zeros(1, nft);
        r = Inf(1, nft);    % initialize residuals

        % traverse each feature in i'th cluster
        for j = 1 : nft
            % feature ft(j) is already in collection?
            idx = ft(j);
            if any(idx == I)
                continue;
            end

            F(:, k + 1) = X(:, idx);
            % [~, ~, t, ~, stats] = regress(Y, F(:, 1 : k + 1), 0.05);
            % r(j)  = sum(t .^ 2) / m;
            % rr(j) = stats(1);
            [b, ~, ~, ~, ~] = regress(Y(1 : s), F(1 : s, 1 : k + 1), 0.05);
            t    = Y(s + 1 : end) - F(s + 1 : end, 1 : k + 1) * b;
            r(i) = mean(t .^ 2);
        end

        if all(r == inf)
            continue;
        end

        % select the feature that has the min residuals
        mini = find(r == min(r));
        mini = mini(1);

        % retraining with complete data
        F(:, k + 1) = X(:, ft(mini));
        [~, ~, t, ~, stats] = regress(Y, F(:, 1 : k + 1), 0.05);

        % udpate collection
        I(k)  = ft(mini);
        R(k)  = mean(sum(t .^ 2));
        RR(k) = stats(1);

        k = k + 1;
    end
end
